Thermodynamics of the two-dimensional Falicov— Kimball model: a classical 

Monte Carlo study 
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The two-dimensional Falicov-Kimball (FK) model is analyzed using Monte Carlo method. In the 
case of concentrations of both itinerant and localized particles equal to 0.5 we determine tempera- 
ture dependence of specific heat, charge density wave susceptibility and density-density correlation 
function. In the weak interaction regime we find a first order transition to the ordered state and 
anomalous temperature dependence of the correlation function. We construct the phase diagram 
of half-filled FK model. Also, the role of next-nearest-neighbor hopping on the phase diagram is 
analyzed. Lastly, we discuss the density of states and the spectral functions for the mobile particles 
in weak and strong interaction regime. 



I. INTRODUCTION 

The study of correlated electron systems has attracted 
great interest over the last decades. Much of this effort 
was devoted to simple Hamiltonians that may contain 
the basic interactions to explain some properties of these 
systems. A common point of departure for many the- 
oretical studies is the Hubbard model 1 which had been 
proposed to describe electron correlations in the narrow- 
band systems. This simple model has been extensively 
investigated over the past forty years, mostly in connec- 
tion with the metal-insulator transition. After discov- 
ery of high temperature superconductors^ it was argued 3 
that the same simple model could possibly capture some 
of the physics of these materials. 

The Hubbard model is represented by the following 
Hamiltonian: 



^Hubb — ^2 t 



nit nil, 



(1) 



where c i(T (c.i a ) creates (annihilates) a conduction elec- 
tron with spin a at lattice site i. The hopping integral 
is usually assumed to be non-zero for nearest neighboring 
sites i and j only. U is the on-site Coulomb interaction. 
Although at first sight such a model may seem oversim- 
plified, only few rigorous results are known. They include 
the one-dimensional solution through the Bethe ansatz, 4 
the Nagaoka's theorem^ and some statements which be- 
come exact in infinity dimensions. 6,7 All other results are 
approximate (mostly of mean-field or perturbative type) 
or obtained for finite lattices, mainly by Lanczos or quan- 
tum Monte Carlo (MC) calculations. 

Therefore, our starting point is a simpler model, that 
can be viewed as a limiting case of a generalized (asym- 
metric) Hubbard model. In the original Hubbard Hamil- 
tonian JQ) the hopping integral t is spin independent. 
However, one can assume that the mass of spin-up and 
spin-down electrons are different, and in the limit of in- 
finitely massive spin-down electrons, they localize and 
only spin-up ones occur in the first sum in Eq. ifl}. Such 
an approximation to the Hamiltonian Q was already 
used by Hubbard i 



Denoting cL (c lT ) by cj (c,), n iT by m and n l]r by 
the resulting Hamiltonian reads 
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Here, Wi is equal to or 1, according to whether the 
site i is occupied or unoccupied by a massive particle. 
The Hamiltonian J2J is known as the FK Hamiltonian. 8 
Within the framework of a common interpretation of the 
FK Hamiltonian, there are two species of particles: itin- 
erant electrons and classical localized particles. The clas- 
sical particles have various physical interpretations: lo- 
calized (/) electrons, spin-down electrons, ions, impuri- 
ties, nucleons. In the following we refer to them as "ions" . 
The ions interact on-site with electrons. There are no 
direct interactions neither between the electrons nor be- 
tween the ions. However, the electron-ion Coulomb in- 
teraction leads to an effective interaction between ions. 
As a result, for a given number of ions, the ground state 
energy depends on their distribution. 

The FK model has a long and successful history in 
dealing with correlated electron systems. Introduced in 
1969 to describe the metal-semiconductor transition in 
SmB6 and related materials? 8 " has been also studied as a 
model of crystallization due to effective interactions me- 
diated by band electronsjSiifi as a binary alloy model and 
many others. The FK model is also useful for describing 
systems that exhibit a phase separatioriii*i2ii&i2ii£ and 
stripe formationiiSiiLi 8 . 

The FK Hamiltonian is over thirty five years old, or 
even older if one takes into account that Hubbard used 
it in 1963 as an approximation to his model. However, 
while it is simpler than the Hubbard model, the general 
solution is also not known. On the other hand, there is 
much more rigorous results for the FK model, then for 
the Hubbard one. One of the most important, proved by 
Kennedy and Lieb^*iS says that at low enough temper- 
ature the half-filled Falicov-Kimball model possesses a 
long range order, i.e., the ions form a checkerboard pat- 
tern, the same as in the ground state. It is a phase, where 
the lattice can be divided into two interpenetrating sub- 
lattices A and B in such a way, that all nearest neighbors 
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of a site from sublattice A belongs to sublattice B and 
vice versa and ions occupy only one of them. This result 
holds for arbitrary bipartite lattices in dimensions d > 2 
and for all values of the interaction strength U. Apart 
from this result not very much is known about solutions 
for two-dimensional systemsii 6 *^*^^! Most of other results 
for the FK model concern one-dimension caseSliSSi2iS4 
or the infinite-dimensional limit The present inter- 
est in this model was renewed in connection with devel- 
oping new calculation methods that can be applied to 
the FK model. In particular, the dynamical mean-field 
theory 2 ^, (DMFT) provides the exact solution of the FK 
model in infinite dimensions. 25 26 While it is known that 
the DMFT captures many key features of the FK model 
even in finite dimensions, this approach also has some 
limitations. The local approximation used in the DMFT 
does not allow to incorporate any nonlocal correlations, 
which are necessary to describe many of phases that are 
expected to realize in the FK model. An extension to the 
DMFT that include short-range dynamical correlations 
has been recently proposedf 2 ^ 2 - 9 Within this approach, 
called dynamical cluster approximation (DC A), the "im- 
purity" used in the DMFT is replaced by a finite-size 
cluster. 

At non-zero temperature the partition function has 
to be calculated for each ionic configuration, and then 
one has to average over the configurations ( "annealing" ) . 
However, since the number of ionic configurations in- 
creases very rapidly with the size of the system, such a 
procedure is possible in some cases only. One of them is 
called "restricted phase diagram" method, 19 where only 
periodic configurations on infinite lattice are taken into 
account. In order to include all possible configurations, 
one has to significantly restrict the size of the lattice. In 
such a case an exact diagonalization technique can be 
appliedf^S This approach is particularly useful for low- 
dimensional systems, where the finite size effects are of 
relatively small importance. Moreover, for linear chains 
it is possible to perform N — + oo extrapolation in a sys- 
tematic way. On the other hand, it is difficult to perform 
the exact diagonalization studies at finite temperatures. 
In order to calculate the partition function one has to 
run a summation over 2 N ionic configurations. There- 
fore, the maximum size of the clusters suitable for exact 
diagonalization study is strongly limited. In Ref. y0| the 
largest clusters consisted of 24 lattice sites. 

In the present contribution we propose to avoid these 
limitations by using the MC method, where thermo- 
dynamic quantities for the FK model arc determined 
by sampling the ionic configuration space stochastically. 
This approach allows us to investigate clusters of up to a 
few hundred sites. Since the frozen particles are classical 
ones, we do not have to use the quantum MC algorithm 
with the "fermionic sign" problem, and thus the calcula- 
tions are not restricted to the high-temperature regime. 

The present calculations are restricted to the neutral 
half-filled case, where the number of ions is equal to the 
number of electrons and their sum is equal to the number 



of lattice sites. The developed formalism can be straight- 
forwardly used away from half-filling, however, the com- 
putational effort in this case is much larger, mainly due 
to the richness of the zero-tempcraturc phase diagram 
of the FK model. Nevertheless, some results for various 
concentrations have already been obtained^ 

The present simulations are performed on a square lat- 
tice. However, also in this aspect the generalization to 
other geometries is straightforward^ 

The outline of the paper is as follows. Section II 
describes the formalism. In particular, it is demon- 
strated there how the classical Metropolis algorithm can 
be adopted to a system with both classical and quantum 
particles. In Section III, we analyze the order-disorder 
phase transition, especi ally in the small-C/ regime. It 
was suggested in Ref. 33 that in this regime the FK 
model exhibits a first order phase transition. Here, we 
analyze this possibility in detail. Section IV is devoted 
to the dependence of the phase-transition-temperature 
on the interaction U. Using temperature dependence of 
the specific heat and charge density wave (CDW) suscep- 
tibility wc identify the critical temperature, constructing 
a phase diagram in the T — U plane. In Section V, we 
demonstrate how the long range order vanishes when the 
temperature increases. Its temperature dependence is 
anomalous in a small-C/ limit. Section VI is devoted 
to the ground state and the phase diagram of the FK 
model with next-nearest-neighbor hopping. Section VII 
presents the results for the density of states and spectral 
functions of the mobile particles. Section VIII contains 
our conclusions. 



II. COMPUTATIONAL METHOD 

In all simulations we have used the Metropolis 
algorithm^ Our system contains classical (ions) as well 
as fermionic (electrons) degrees of freedom. The appro- 
priate way to treat such a Hamiltonian is to define the 
grand canonical partition function as 

Z = ^Tr e e-^ c )-H, (3) 
c 

where TC(C) is the Hamiltonian @ for a fixed ionic con- 
figuration C, indicates summation over ionic config- 
urations, Tr e denotes the trace over fermionic degrees of 
freedom, /3 is the inverse temperature, and N is the oper- 
ator for the total number of electrons. For a given ionic 
configuration the Hamiltonian Ti(C) can be diagonalizcd 
numerically and the summation over fermionic degrees 
of freedom gives the partition function in the following 
form 

Z = Y,H{l + e-e [B ^ e) -» ] }, (4) 

C n 
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where E n (C) are eigenvalues of H(C). Introducing the 
electronic free energy 

Fe(C) = -~J2ln{l + e-^ c ^}, (5) 

n 

the partition function can be written in a form analogous 
to that used for a spin system 

c 

The above equation describes the effective model, that 
results from the summation over the fermionic degrees of 
freedom. As the electronic free energy T e {C) depends on 
temperature, the corresponding Hamiltonian, describing 
only the classical particles would include temperature- 
dependent interaction. This dependence may lead to 
nontrivial critical exponents. Equation © indicates also, 
that in the Metropolis algorithm we should use the elec- 
tronic free energy in the statistical weights. To do this, 
however, one has to know the value of the chemical po- 
tential /i, and apart from specific cases (e.g., half-filling) 
it has to be determined separately. Since the calcula- 
tions are carried out in the grand canonical ensemble, 
the chemical potential must be kept constant during the 
simulation. Determining thermodynamic quantities, the 
averages are calculated with the statistical weights 

w{C) = ±e~^ c \ (7) 

of corresponding ionic configurations C. 

One of the advantages of the proposed approach is 
that it gives densities of states and spectral functions 
for the mobile particles that identically satisfy different 
sum rules. This will be discussed in Sec. VII. 

The calculations were carried out on clusters up to 
10 4 lattice sites, however, in most cases square clusters 
20 x 20 were used. Usually, the simulations started at 
high temperature, with the initial state with randomly 
distributed ions. Then, the temperature has been slowly 
decreased. Below some temperature, ions started to form 
a checkerboard patterns. 

It is difficult to describe the process of pattern forma- 
tion qualitatively. In some cases it is convenient to use 
A = \{ua) — («b)|) where (ua) (("b)) is the ion con- 
centration in sublattice A (B), as the order parameter. 
However, such a description is not well suitable for cluster 
calculations, since even in the case, when in the ground 
state ions form almost perfect checkerboard, A can be 
close to zero. Figure 1 presents an example of such a 
ionic configuration. A small staggered field can be intro- 
duced in order to prevent an occurrence of such effects. 
However, this "phase smoothing out" may restrict also 
some other configurations. 

Instead, we use the density-density correlation func- 
tion for the ions to describe the ordered state. It has 
the advantage that it is capable to describe long range as 



FIG. 1: Example of almost perfect checkerboard ionic config- 
uration, for which the long range order parameter A is close 
to zero. Filled circles represent sites occupied by the massive 
particles. 

well as short range correlations. We define the correlation 
function in the following way: 

I N 

5n= 4A?H Yl w{r i )w{r i +T 1 x + T 2 y), (8) 

i— 1 ri,T2=±n 

where w (rj) = Wi, x, and y denotes unitary vectors along 
the x and y directions, respectively (the lattice constant 
a = 1 has been assumed). Note, that g n describes corre- 
lations along the lattice axes. 

The correlation function g n describes ionic configura- 
tion only. However, the distribution of the light particles 
correlates with the ionic distribution due to the direct in- 
teraction between these two types of particles. The larger 
value of U the stronger correlations take place. Figure 2 
demonstrates these correlations for weak (upper config- 
urations) and strong (lower configurations) interaction. 

If there is a phase transition, some thermodynamical 
quantities have to diverge at the critical temperature. 
Here, we used the CDW susceptibility x an d the specific 
heat Cy to determine this point. The CDW susceptibility 
is related through the fluctuation-dissipation theorem to 
the variance of the density-density correlation function. 
This is a form especially convenient for the MC calcula- 
tions 

x = ^ f (( 9 2 1 )-(gi) 2 ), (9) 

where we used (. . .) to indicate the average over gener- 
ated ionic configurations. It is a little bit more compli- 
cated to determine the specific heat. Within the stan- 
dard MC approach the specific heat is calculated from 
the fluctuations of the energy, similarly to Eq. JjJ}. In 
our case, however, according to Eq. (J7J) we use the elec- 
tronic free energy in the Metropolis algorithm and the 
internal energy is not directly available from the simula- 
tions. On the other hand, the Fermi energy is much larger 
than the order-disorder transition temperature. There- 
fore, estimating the specific heat we replace the trace over 
the fermionic degrees of freedom by a ground state ex- 
pectation value E, 35 which is a temperature-independent 
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FIG. 3: Specific heat determined from Eq. 1101 (C v , dashed 
line) and from Eq. (C|°, solid line) for U/t = 1. For 
comparison, CDW susceptibility \ is a l so presented (dashed- 
dotted line). 



FIG. 2: Examples of ionic (left) and electronic (right) config- 
urations for weak (U/t — 0.5, fcgT/t = 0.1, upper row) and 
strong (U/t = 5, k-gT/t = 0.3, lower row) interaction. Diam- 
eters of the circles in the right configurations are proportional 
to the concentration of electrons. 



quantity. Then, the specific heat can be determined from 
the fluctuation-dissipation theorem 



1 1 

N k B T' 2 



((E 2 ) (E) : 



(10) 



In order to confirm the validity of this approxima- 
tion we compare the specific heat determined numerically 
from the relation 



dU 
dT 



(11) 



and that from Eq. (ID) . Since the simulations were car- 
ried out only for a given set of temperatures, a finite 
difference approximation had to be used for the deriva- 
tive in Eq. (|1 1|) . Nevertheless, within the temperature 
regime of interest the difference is below a few percent. 
Figure |3| shows a comparison of the specific heat calcu- 
lated from Eq. lfTT)|) and from Eq. fTTJ) . Additionally, 
the CDW susceptibility is plotted there in order to indi- 
cate the critical temperature. A comparison of the crit- 
ical temperatures determined from the specific heat and 
from the CDW susceptibility, presented in Figs. ijTtjHT^)! . 
confirms the validity of the proposed approach as well. 
Note, that these quantities are calculated in completely 
different ways: the specific heat is determined from the 
eigenvalues of the fermionic Hamiltonian, whereas the 
CDW susceptibility from the correlation function of the 
classical particles. 

In each simulation the ensemble averages of the ther- 
modynamic quantities of interest are calculated after the 
system equilibrated. The thermalization period varies 
with temperature. In particular, it is very long at low 



temperatures or in the vicinity of critical points (the crit- 
ical slowdown). In order to determine whether the sys- 
tem has reached equilibrium, we usually start the simu- 
lation with two replicas of the system: one starting from 
a fully ordered state, the other from a random ion distri- 
bution. Then, the energies of the systems evolve crossing 
each other after some number of MC steps. This point 
is considered as the beginning of equilibrium state and 
the averages of the thermodynamic quantities are calcu- 
lated over the remaining MC steps. This procedure is 
especially useful close to the transition temperature. 

The width of the peak in the specific heat and in the 
CDW susceptibility decreases with the increase of the 
size of the system and a finite-size scaling should be 
carried out in order to determine precisely the transi- 
tion temperature. This can be done, for instance, using 
the standard Binder cumulant method. 36 However, tak- 
ing into account the huge amount of time needed to run 
the finite-size scaling over the whole parameter space, we 
have decided to omit to do it and most of the presented 
results have been obtained for 20 x 20 cluster. Such a 
cluster is sufficiently large to produce results accurate 
enough to describe the properties of the model under in- 
vestigation and, on the other hand, the time of a single 
simulation run is short enough to allow determination 
of the full phase diagram. Moreover, some comparison 
of results obtained for 20 x 20 and 40 x 40 clusters is 
presented in Section VII. 



III. NATURE OF THE PHASE TRANSITION 

For large U there is an Ising-like phase transition at 
the critical temperature T c oc U . In the small-J7 
regime the critical temperature is bounded from below 
by U 2 /\ lnJ7|i2iiS In fact, it was shown that in large di- 
mensions T c oc U 2 \ In U\ for U -> 

It is known that for large U the Falicov-Kimball model 
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belongs to the same universality class as the Ising model 
and the order parameter is described by Curie- Weiss law 
A = tanh(A/6), where 6 = T/TJJJ). Thus, the phase 
transition is of second order in this regime. On the other 
hand, it was shown in Refs. l26l and l3^ that in infinite di- 
mensions in the small-?/ limit the order parameter has a 
strange non-BCS-like temperature dependence. There- 
fore, it is important to determine precisely the nature 
of the phase transition in this regime. We have used a 
method proposed by M. Challa et al£& to distinguish be- 
tween first and second order phase transitions. Systems 
undergoing first order phase transitions are accompanied 
by free energy barriers which separate the free energy 
minima characterizing the coexisting phases. It results 
in discontinuities in the first derivatives of the free en- 
ergy, e.g., the internal energy. This, in turn, leads to a 
(5-type singularity in the specific heat at the transition. 

Within the framework of Metropolis MC approach 
the energy E fluctuates with the probability distribution 
P{E) usually given by a Gaussian. Its width is propor- 
tional to the specific heat. However, if the system is close 
to the first order transition, the probability distribution 
P(E) is a superposition of two Gaussians centered at 
different energies, E + and Here, E + and £L are 

the energies in the high- and low-temperature phases, 
respectively. The Gaussians are weighted by the Boltz- 
mann factors of E + and and thus this splitting oc- 
curs in a vicinity of the transition temperature only. At 
higher (lower) temperatures P{E) forms a single Gaus- 
sian centered at E + (E-). 

Figure 01 shows the energy distribution for tempera- 
tures slightly below and above the transition temperature 
for U/t — 0.5. Since the data has a double-peak struc- 
ture, the phase transition is of the first order. One can 
see the transfer of the weights from the low-temperature 
Gaussian peak (the lower panel in Fig. 0} to the high- 
temperature one (the upper panel), as the temperature 
increases, passing through the critical value. The energy 
distribution at the critical point, where the heights of 
both the peaks are comparable, has been presented in 
Ref. |H 

At temperatures much lower and much higher than 
the transition temperature, the energy distribution can 
be well fitted by a single Gaussian. Such situations are 
presented in Fig. 

The energy distribution at low temperature (lower 
panel in Fig. 0) consists of a large number of peaks, 
clearly visible in the inset. They are connected with the 
discrete spectrum of the Hamiltonian for a finite system. 
Namely, at such a low temperature, there are only a few 
dislocations in the checkerboard ionic configuration. Due 
the low concentration of the dislocations, they are almost 
independent and each of them changes the energy by an 
approximately the same amount. In this way, two succes- 
sive peaks correspond to configurations with the numbers 
of dislocations that differ by one. As the number of dislo- 
cations increases, they start to "feel" each other, and the 
effective interactions smear out this energy ladder. This 




FIG. 4: Probability distribution of the energy at temperatures 
close to the critical temperature for U/t = 0.5. The thick 
solid lines represent a superposition of two Gaussians that 
fit the simulation results, whereas the dashed and dashed- 
dotted lines show the component Gaussians. The arrows in- 
dicate positions of the centers of the Gaussians representing 
low-temperature (E-), and high-temperature (E+) phases, 
respectively. 



is why the peaks are visible only in the low-energy part 
of the distribution. Of course, when the size of the lattice 
increases, the energy spectrum becomes quasicontinuous 
and the oscillations disappear. 

The coexistence of low- and high-temperature phases 
at the phase transition can also be observed in the CDW 
susceptibility. Fig. shows the two-peak structure of 
the probability distribution of x- It should be noted that 
the distribution presented in Fig. describes ionic con- 
figurations, whereas the one presented in Fig. 0] is ob- 
tained from the eigenvalues of the fermionic Hamiltonian. 
The similarity between these two distributions speaks 
strongly in favor of the validity of the proposed MC ap- 
proach, confirming the presence of first order phase tran- 
sition in the small-t/ regime. 

The magnitude of splitting of the energy distribution 
close to T c decreases with the increase of the interaction 
strength U and above a critical value U* disappears. The 
same holds true for the distribution of the CDW suscep- 
tibility. Figure presents the energy and CDW suscep- 
tibility distribution at the phase transition for U/t = 3. 
One can notice an excellent agreement with the theoret- 
ical Gaussian curve, indicating an absence of any phase 
coexistence and the second order character of the phase 
transition. 
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FIG. 5: Probability distribution of the energy at temperatures 
much higher and much lower than the transition temperature 
for U/t = 0.5. Inset in the lower panel shows the low-energy 
part of the distribution with lines connecting the points. 
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FIG. 6: Probability distribution of the CDW susceptibility \ 
close to the phase transition for U/t = 0.5. The meaning of 
the lines is analogous to that of Fig. [I] 

We have estimated the critical value [/*«t, however, 
extensive numerical studies are necessary in order to de- 
termine U* precisely. 



A. Phase diagram 

The position of peaks in the specific heat and the CDW 
susceptibility has been used to determine the critical tem- 



0.10 




E/t X 



FIG. 7: Probability distribution of the energy (left panel) and 
of the CDW susceptibility (right panel) at the phase transition 
for U/t = 3. 

perature. Fig. [21 shows the typical temperature depen- 
dence of the specific heat. The phase diagram has been 
constructed using results obtained for 20 x 20 system. 
Figure [8] presents the transition temperature as a func- 
tion of the interaction strength. The horizontal axis is 
plotted as U/(U + t), such a scaling allows one to present 
both the weak and strong coupling results in the same 
graph. There are also shown results obtained by means 
of other methods: DMFT for d = oo case, taken from 
Ref. Ejl MC taken from Ref. |U and DCA taken from 
Ref. |23- In the large-£7 limit the critical temperature for 
the Ising model is also presented for comparison. 



IV. ORDER PARAMETER 

In order to quantitatively describe the ionic configura- 
tion, we have investigated the density-density correlation 
function defined in Eq. (|SJ). This function is capable of 
descring short range as well as long range order. In a 
fully ordered (checkerboard) state it oscillates with a pe- 
riod equal to two lattice constants (see Fig. |5J). 

It is convenient to define a renormalized correlation 
function G n = (— l) n 4 (g n — pf), pi = N. t /N is the con- 
centration of ions, where N and 2V$ are numbers of lattice 
sites and ions, respectively. Such a function is equal to 1 
for the checkerboard state and close to for randomly dis- 
tributed ions, independently of the distance n (for n > 0). 
Apart from these limiting cases, this function decreases 
monotonically with increasing distance. Figures HOtfO 
present the temperature dependence of G n for a wide 
range of the interaction U . These curves are presented 
together with the temperature dependencies of the spe- 
cific heat and the CDW susceptibility, calculated for the 
same values of U. This allows one to find the exact values 
of the critical temperatures. 

For strong to intermediate values of the Coulomb in- 
teraction there is a distinct peak in the specific heat, 
indicating the phase transition from the ordered state to 
the disordered one. The corresponding vanishing of the 



7 




0.4 0.6 

U/(U+t) 

FIG. 8: Transition temperature for the d = 2 half- 
filled Falicov-Kimball model as a function of the interaction 
strength (solid line with circles) . The dotted and dashed lines 
present for comparison results taken from Ref. W\ These 
results were obtained in d = oo limit for the Bethe lattice 
(dotted line) and the hypercubic lattice (dashed line). Open 
circles represent numerical results taken from Ref. l2fiL esti- 
mated as the temperature at which a gap in the density of 
states closes. Stars represent results obtained in Ref. [2^ from 
DCA with a QMC for a 36-cite cluster. The solid line con- 
necting filled circles is a guide for the eyes only. 




kj/t 



FIG. 10: Specific heat Gy and CDW susceptibility \ (upper 
panel) and correlation function G n (lower panel) as a func- 
tion of temperature for U/t = 20. Various lines in the lower 
panel correspond to correlation functions G n calculated for 
various distances n. For comparison, there is also a line rep- 
resenting temperature dependence of the magnetization in the 
Ising model 
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FIG. 9: Correlation functions g n (dashed lines) and G n (solid 
lines) for U/t = 20. Left panel shows the case of a fully 
ordered (checkerboard) state and the right one corresponds 
to the vicinity of the phase transition. 



correlation function G n resembles a typical behavior of 
an order parameter close to the second order phase tran- 
sition (see Figs. 1101 and ITD) . Note, that the long range 
correlations disappear more rapidly than that for shorter 
distances. This can be, however, attributed to the finite 
size of the system. 

On the other hand, in the weak interaction limit 
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FIG. 11: The same as in Fig. Holbut for U/t = 1. 
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FIG. 12: The same as in Fig. ITUlbut for U/t = 0.01. Addi- 
tionally, there is also a line representing the order parameter 
determined from DMFT solution of the FK model in infinite 
dimensions in the limit U — » 0. This solution is taken from 
Ref. HI 



the correlation functions vanishes almost linearly and 
the slope begins already at temperature approximately 
0.25T C . An unusual behavior in this regime is seen also in 
the specific heat: there is hump in both Cy(T) and x(T) 
close to the temperature at which the slope in G n (T) oc- 
curs. Fig. 1121 illustrates such a behavior for U/t = 0.01. 
For weaker interaction the values of the specific heat and 
the CDW susceptibility are different, however, the shape 
of their temperature dependence, as well as the temper- 
ature dependence of the correlation function, remain al- 
most unchanged. 

Such an anomalous temperature evolution of G n (T) 
may be explained as a result of finite size of the system. 
For small values of U there exist solutions with periodic 
ionic configurations possessing very large period. Such 
configurations are excluded from our calculations due to 
limited size of the clusters the simulations were carried 
out on. 

On the other hand, it was recently shown that similar 
behavior of the order parameter occurs in the infinite- 
dimensional limiti2L22i Our results may suggest that such 
a behavior is a generic property of the weak-coupling 
Falicov-Kimball model. 
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FIG. 13: Possible ground-state ionic configurations in t' 2> t 
limit. The dashed and dotted lines indicate (independent in 
t — > limit) sublattices. 



V. NEXT-NEAREST-NEIGHBOR HOPPING 

A. Ground state 

In this section we generalize the FK Hamiltonian by 
taking into account the hopping to next nearest neigh- 
bors (NNN) with the hopping integral t' in addition to 
the nearest-neighbor (NN) hopping. In particular, the 
question concerning the nature of the ordered state is 
addressed. In a more general case of various concentra- 
tion of both type of particles this problem was analyzed 
in Ref. E 

In the limiting case of t' <C t one may expect that 
the checkerboard state is still the actual ground state of 
the FK model. In the opposite limit, when the hopping 
to the nearest sites can be neglected, the square lattice 
can be divided into two interpenetrating square sublat- 
tices with the lattice constant V2 times larger than the 
original one and with the axes rotated by 45°. Since 
the electrons do not hop between the sublattices, the 
system breaks up into two uncoupled, interpenetrating 
square lattices composed of NNN bonds. Each of these 
lattices is independently described by the FK Hamilto- 
nian. Since both the sublattices are bipartite and both 
the subsystems are half-filled, the ions will arrange into 
checkerboard patterns in the ground state. Depending 
on the relative phase between the orderings in the sub- 
lattices, the resulting ground state of the whole system 
will have the form of vertical or horizontal stripes. Fig. 
1131 demonstrates these possibilities. In the large-?/ limit 
the threshold value of the ratio t'/t, that separates the 
ground states with ions forming the checkerboard and 
stripe patterns, can be determined using the mapping of 
the FK model onto the Ising one. In the antiferromag- 
netic Ising model with NN and NNN coupling the ground 
state is the simple antiferromagnet for J' / J < 0.5, where 
J and J' are the NN and NNN interactions, respec- 
tively. For J'/ J > 0.5 the system minimizes the energy 
by ordering in alternate ferromagnetic rows of opposite 
spinsi 41 i 42 i 43 Such spin configurations ( "superantiferro- 
magnetic" ) correspond to the ionic configurations of the 
FK model presented in Fig. Since the ratio J'/ J of 
Ising NN and NNN interactions is equal to (t'/t) 2 in the 
corresponding FK model, the threshold value of the NNN 
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FIG. 14: Specific heat for U/t = 1 and t'/t = 0, 0.3 and 0.5. 

hopping is given by t'/t ~ 0.71. In Ref. |40| it has been 
shown that this result holds true up to fourth order in 
t/U. For weaker interaction this threshold ratio can be 
determined from a comparison of the ground state energy 
of both the checkerboard and stripe configurations. Our 
simulations indicate that the critical value oft'/t changes 
slightly with U , however, it always lies between 0.71 and 
0.8. 



B. Specific heat 

Since the NNN hopping term introduces frustration, 
it reduces the temperature of the transition from the 
checkerboard to the disordered state (for t'/t < 0.71). 
In the strong interaction limit the dependence of the 
critical temperature on t'/t can be determined from re- 
sults obtained for the corresponding Ising model with 
NNN interaction. For finite U/t the critical temperature 
can be identified from the position of the peak in the 
temperature dependence of the specific heat. Figure PHI 
presents results for U/t = 1. The critical temperature 
decreases with increasing t', going to zero for t'/t ~ 0.71. 
At this point there is no long range order at any finite 
temperature. Then, the critical temperature increases 
with further increasing t' . In this regime the peak in 
the specific heat indicates the transition from the stripe 
configuration to the disordered state. For t'/t ^> 1 the 
nearest-neighbor hopping can be neglected and one ends 
up with the FK model on a square lattice with NN hop- 
ping integral t'. Therefore, according to the phase dia- 
gram presented in Fig. [SJ in the limit t'/t — > oo (for a 
given U) the critical temperature goes to zero. Figure IT5l 
shows how the specific heat depends on t'/t at tempera- 
ture k^T/t = 0.02. This is a critical temperature for two 
different values of t'\ one for t'/t = 0.4, when in the low- 
temperature phase the ions form the checkerboard pat- 
tern, and the other for t'/t = 1.2, when stripes minimize 
the free energy at low temperature. As a result, the spe- 
cific heat plotted as a function of t'/t has two maxima, as 
can be seen in Fig. El The left peak corresponds to the 
transition from the checkerboard state, whereas the right 
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FIG. 15: Specific heat for U/t = 1 as a function of t'/t at 
k^T/t — 0.02. The arrows indicate values of t' , for which 
ksT/t — 0.02 is the critical temperature of the phase tran- 
sition from the checkerboard state (A) and from the stripe 
phase (B). The inset shows temperature dependence of the 
specific heat for t'/t — 0.3 and t'/t = 1.4. 

one to the transition from the stripe state. The widths of 
the peaks are connected with the effective sizes of the lat- 
tices. For large NNN hopping the lattice breaks up into 
two sublattices, which becomes decoupled in the limit 
of t'/t — > oo and the checkerboard pattern is formed in 
each of them independently. Therefore, the simulations 
are effectively carried out for two replicas of a system of 
halved size. As a result, the system is further away from 
the thermodynamic limit than in the case of small t' and 
the maximum of the specific heat is less sharp (see the 
inset in Fig. IT5|) . The features seen at the top of the sec- 
ond maximum are connected with a rearrangement of the 
ions from two-dimensional structures to one-dimensional 
ones, as can be seen in the lowest row of Fig. El Figure 
1161 illustrates how the ionic configuration evolves from 
the checkerboard to stripe pattern with increasing am- 
plitude of the NNN hopping. Due to the degeneracy of 
the ground state for t'/t > 0.71, showed in Figure El 
the evolution goes with equal probability through a state 
with vertical or horizontal stripes. Finally, for any finite 
temperature the system ends up in a disordered state. 



C. Phase diagram 

The determination of the full phase diagram of the FK 
model with NNN hoppings is a very CPU time consuming 
task due to the large number of free model parameters. 
Therefore, we have determined only a few points in order 
to qualitatively describe the dependence of the critical 
temperature on the ratio t'/t. Figure El shows the re- 
sults for U/t = 1. This figure does not show the decrease 
of the critical temperature, when t' becomes very large. 
This is because we show the temperature in units of the 
NN hopping integral t. In the case of small NNN hopping, 
t is directly connected to the band width. However, when 
t' is very large, the band width is connected to t' rather 
then to t and k-gTc/t', instead of fceTc/i, would go to 
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FIG. 16: Snapshots of ionic configurations for various NNN 
hopping amplitudes. The simulations were performed for 
U/t = 1 and k B T/t = 0.02. 
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FIG. 17: Critical temperature as a function of t'/t for U/t = 
1. The line connecting points is a guide for the eyes. The 
configurations indicate the ground states in given t'/t regimes. 
The inset shows the critical temperature in U/t — ► oo limit. 



zero when t' goes to infinity. The inset demonstrates the 
corresponding phase diagram in U/t — ► oo limit obtained 
through mapping onto the NNN Ising model. Since the 
critical temperature goes to zero when U goes to infin- 
ity (for a given bandwidth), the temperature has been 
plotted as k B TU/t 2 . 



VI. DENSITY OF STATES 

Since the electron distribution correlates with the ionic 
configuration, the formation of the checkerboard pattern 
results in a modification of the electronic density of states 
(DOS). The DOS at a given temperature can be obtained 
by averaging densities determined in each MC sweep. For 
weak electron-ion interaction at high temperature the 
electrons on average are almost unaffected by the ionic 
configurations and its DOS resembles that for free elec- 
trons on square lattice, except for the lack of the van 
Hove singularity (see Fig. 118b). As the temperature is 
lowered the ions form the checkerboard pattern and the 
electronic concentration decreases in sites that are oc- 
cupied by the ions. As a result of this CDW order a 
gap opens m the DOS at the Fermi level (see Figs. HHk 
and 118b ). The corresponding metal-insulator transition 
remains in accordance with the Mott pictured The sit- 
uation is different in the strong interaction regime. In 
this case, the electron-ion interaction almost completely 
forbids electrons from occupying sites that are already 
occupied by the ions. Therefore, even at high tempera- 
ture the DOS averaged over the MC sweeps differs from 
the free electron DOS. At high enough temperature the 
Metropolis algorithm accepts every ionic configurations 
with probability almost equal to one. In other words, 
at high temperatures the ionic configurations are not af- 
fected by the electron-ion interaction and the ions are 
distributed randomly. In this regime, the FK model re- 
duces itself to a model of free electrons in a random po- 
tential of disordered binary alloy. For large enough U 
there is a gap in the DOS, that does not change further 
with the increasing temperature (see Fig. 119b). When 
the temperature is lowered, this gap transforms itself into 
the CDW gap (Figs, and El)). 

There is a significant difference in the temperature evo- 
lution of the DOS between the FK model and a CDW 
system with selfconsistently determined gap»^ In the 
FK model the distance between the edges of the gap is 
constant when the temperature increases, but its depth 
as well as the high of the peaks decreases. It will be 
clearly visible in the proceeding section, where the spec- 
tral functions are presented. On the contrary, the width 
of the selfconsistently determined CDW gap decreases 
with increasing temperature, whereas its depth remains 
constant. 

The DOS presented in Figs. 1181 and IT§1 are obtained 
for 20x20 clusters. The question arises whether the fine 
structure seen in the DOS is a finite-size effect or is in- 
herent to the FK model. The finite-size effects are es- 
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FIG. 18: Density of states for U/t = 1 at various temperatures. For temperatures T higher then fcgT/t = 0.1 the DOS is 
almost the same as in panel c). The inset in panel c) shows the corresponding DOS for 40 x 40 cluster. 
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FIG. 19: Density of states for U/t = 8 at various temperatures. For temperatures T higher then k-gT/t — 1 the DOS is almost 
the same as in panel c). 



pecially important for small U. This drawback is visible 
in Fig. 118b . where the inset shows the DOS obtained for 
40x40 cluster. In both cases of 20x20 and 40x40 clus- 
ters, the overall DOS resembles corresponding DOS for 
free electrons on square lattice. However, for the DOS 
obtained for 40x40 cluster is much more smooth, and 
therefore we attribute the fine structure to the finite size 
of the cluster. The DOS presented in Fig. 118b does not 
change with increasing temperature. This indicates that 
for U/t = 1 temperature kBT/t — 0.1 is high enough 
to treat the electrons described by the FK Hamiltonian 
as a free electron gas on a square lattice with diagonal 
disorder. Moreover, the interaction U is weak enough to 



neglect in this temperature regime the influence of aver- 
aged disorder and the actual DOS is almost the same as 
for free electrons on a square lattice. 

On the other hand, when U increases, the DOS be- 
comes much less dependent on the cluster size and the 
features that develop in the DOS can be attributed to 
the FK model itself. In order to confirm this assump- 
tion we compare DOS calculated for 20x20 cluster with 
delta function broadening r] — 0.02 1 and for 40x40 with 
r\ = 0.005 t. Results are presented in Fig. |201 The 
smaller broadening used in the case of simulations on the 
larger cluster can uncover more details of the actual den- 
sity of states. The differences between the DOS obtained 
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FIG. 20: Comparison of DOS obtained for 20x20 and 40x40 
clusters with different delta function broadening. Note, that 
the right vertical axis is shifted in order to make the figure 
more legible. The presented results have been obtained for 
U/t — 8 at temperature ksT/t — 1, i.e., the same as in Fig. 

nit. 



for 20x20 and 40x40 clusters with the same broadening 
7] = 0.02i are within the linewidth. The results presented 
in Fig. |20l can be directly compared to results obtained 
using dynamical cluster approximation (see Fig. 6 in Ref. 
128^ . When the size of the cluster increases from lxl 
(what corresponds to the dynamical mean-field approx- 
imation) to 8x8, some features start to develop in the 
DOS. The positions of these features are almost the same 
as in our MC calculations, however, they are less devel- 
oped. In particular, it seems that the DOS at u> = ±U/2 
vanishes in infinite system and narrow peaks are located 
at \lu\ slightly larger than U/2. In Ref. |2|| these peaks 
are attributed to the localization of electrons in four sites 
surrounding each site occupied by an ion. However, these 
features, as well as the overall density of states, do not 
change when the temperature increases. Moreover, the 
high-temperature density of states is exactly the same 
as the one for systems with random on-site binary disor- 
der. Both the peak and the (pseudo)gap in the vicinity of 
±E//2 have been obtained within many approaches to dis- 
ordered binary alloys4& This indicates that those peaks 
are rather connected with the strong interaction between 
electrons and randomly distributed ions. 



VII. SPECTRAL PROPERTIES 

The temperature and interaction dependence of the 
density of states is connected with the spectral proper- 
ties of the FK model. Since at half filling there is no 
phase separation, the system is translation invariant and 
the ionic-configuration-averaged electronic Green func- 



J2 ex P ^ {k-Ri-kf- Rj)} (G {Ri,R h z)) 

= G(k,z)5(k-k'). (12) 



Ri Rj 



Here, Q (Ri,Rj, z) — <[z — H(C)] > is the real-space 
I J ij 

Green function for a given ionic configuration C and (...) 
denotes averaging over the configurations. Then, the 
spectral function can be determined from the standard 
formula 



A(k,w) 



1 



ImG(fc,cj + £0 + ) 



(13) 



Figure [2] shows the spectral weights corresponding to 
the DOS presented in Figures ^] and ^| In the rela- 
tively weak interaction regime (U/t — 1) at low tem- 
perature there is the CDW gap at (0, n) and (tt/2,tt/2) 
(Fig. 121b). The dispersion relation is very accurately 



described by — ±y (efe — /f) + U 2 , where tk = 
— t(cosk x + cos k y ) is the energy of free electrons on a 
square lattice. With the increasing temperature the ionic 
disorder increases, destroying the checkerboard pattern. 
As a result, the spectral peaks gradually broadens up fill- 
ing the gap (Fig. 121b ). However, due to the weakness 
of the electron-ion interaction, this broadening is rela- 
tively small. Eventually, at high temperature the CDW 
gap disappears and the dispersion relation takes on the 
form of efc, that does not change with further increase of 
temperature (Fig. 121b). In the strong interaction regime 
(U/t = 8) at low temperature the dispersion relation can 
be described by the same formula, as for U/t =1. In this 
case, however, the stronger interaction leads to almost 
equal spectral weights in the lower and upper subband 
and the band structure does not resemble that for free 
electrons any more (Fig. I21t0. When the ionic disorder 
increases, additional bands appear, mainly around (0, 0) 
and (it, it) points (Fig. 121b). This effect is more visible in 
the corresponding DOS (Fig. IT§b). These bands broaden 
up with farther increase of temperature and most of the 
spectral weight is transferred to them (Fig. 12 It). Despite 
the broadening, at point (0, 7r) there still exist flat parts 
of the bands, that are responsible for the peaks in the 
DOS at energy ±U/2. 

There exist some exact results for the mobile particles 
in the form of sum rules and the numerical results can 
be checked against them4L22i They give values of a few 
lowest spectral moments, defined by 



wM(k,w) 



(14) 



Their values in the half-filling case are given by very 
simple expressions 



jU (k) = 0, /ii(k) = e k , /x 2 (k) = e 2 . 



U 2 



(15) 
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FIG. 21: (Color online) Spectral functions for U/t = 1 (upper row) and U/t — 8 (lower row) at various temperatures. The 
dotted line indicates dispersion of free electrons on the square lattice. 



It is interesting that the same results are valid for the 
Hubbard models and for the FK model in a nonequi- 
librium stated Since in the proposed approach the FK 
Hamiltonian for a finite system is exactly diagonalized in 
each MC step, all these sum rules are exactly satisfied. 
As a result the sum rules are also satisfied for the spectral 
functions obtained in the whole MC run. However, the 
results for the moments are exact only if they are deter- 
mined directly from the distribution of the eigenenergies 
in each MC step. The other approach, i.e., integrating 
u! n with the spectral function obtained in the whole MC 
run, leads to a small error introduced by the broadening 
of the Dirac 8 functions. The same holds true for the 
local moments^ where w™ is integrated with the density 
of states. 

VIII. SUMMARY 

We have presented Monte Carlo analysis of half-filled 
FK model. In order to take into account both the classi- 
cal and fermionic degrees of freedom, we have derived a 



modification of the classical Metropolis algorithm, where 
the interaction energy is replaced by free energy, calcu- 
lated by numerical diagonalization of the FK Hamilto- 
nian for a given ionic configuration. Such an approach is 
possible due to the absence of many-body interactions in 
the FK Hamiltonian. Although there is no many-body 
term in the Hamiltonian, averaging over ionic configu- 
rations ("annealing") leads to many-body effects in the 
FK model. As a result, the FK model possesses a rich 
phase diagram. The simulations presented in this paper 
concern the case of pi = p e = 0.5, i.e., the case where the 
ground state is known to be the checkerboard ionic con- 
figuration. Our results illustrate how the system reaches 
this state when the temperature is lowered. It is known, 
that in the strong interaction limit the FK model maps 
onto the Ising model, and therefore the phase transition 
to the ordered state is of second order. There is no such a 
rigorous result for weakly interacting 2D FK model. Our 
simulations indicate the presence of a first-order phase 
transition in this regime. However, the precise determi- 
nation of the critical value of the interaction, below which 
the phase transition is of first order, requires extensive 
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simulations, mostly in order to perform the finite-size 
scaling. 

Additionally, we have shown that the order parame- 
ter decreases unusually with increasing temperature in 
this regime. A departure from the Ising-type tempera- 
ture dependence of the order parameter has recently been 
demonstrated also in infinite dimensions^ 

Performing the MC simulations for various interac- 
tion strengths and at various temperatures, we have con- 
structed the phase diagram for half-filled FK model. It 
was shown, that the boundary line, separating the or- 
dered (checkerboard) and disordered phases, determined 
from positions of the peaks in the specific heat and in 
the CDW susceptibility coincides with the one obtained 
from the opening of a gap in the electronic DOS. The 
phase diagram has also been determined in the presence 
of NNN hopping. It was shown, that there is a point (for 
t'/t w 0.7 -j- 0.8, depending on the interaction strength) 
at which no long range order exists at any finite tem- 
perature. For larger t' the ions at low temperate form 
horizontal or vertical stripes. 

The electronic DOS and spectral functions have been 
determined for a wide range of the interaction strength. 
In the large-?/ limit, the distribution of electrons strongly 
depends on the ionic configuration. As a result, the CDW 
gap develops in the DOS at low temperatures. At high 
temperatures, the ions are distributed randomly over the 
lattice sites and the electronic DOS is the same as in a 
free electron system with a diagonal disorder. On the 
other hand, in the weak interaction regime the electrons 
are hardly affected by the ions, and averaging over the 



ionic configurations leads to the same band structure as 
for free electrons on a square lattice. We have shown, 
that the proposed approach give both the electronic DOS 
and spectral functions in a good agreement with that ob- 
tained with use of much more elaborated DCA method, 
whereas the computational effort is much lower. 

The present paper demonstrated how the MC ap- 
proach can be used to describe the formation of the 
checkerboard order in the half-filled FK model. It is also 
possible to use this method to analyze the FK model 
away from half-filling, where inhomogeneous states are 
expected to minimize the free energy at low enough tem- 
peratures. The MC method would be particularly use- 
ful to describe phases with irregular ionic configurations, 
e.g., phase separation, which are difficult to analyze ana- 
lytically. Preliminary result in this area are presented in 
Ref. 31. Another area where the Monte Carlo study of 
the Falicov-Kimball may be particularly useful, includes 
nonbipartite lattices. In this case, the checkerboard pat- 
tern cannot be formed even at half-filling, and the model 
is expected to exhibit a strong frustration. This issue is 
under current investigation. 32 
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